# %%
import pandas as pd
import scanpy as sc

# %% the mtx output from read_mtx() are often need to be transposed to make var = gene, obs = cell
# read_mtx & related func cannot interpret ~ as user dir?
mtx = sc.read_mtx('data/ren2021/matrix.mtx.gz').X.transpose()

adata = sc.AnnData(mtx)
adata

# manually read barcodes & genes
cell = pd.read_table('barcodes.tsv.gz', names = ['id'])
cell

adata.obs_names = cell['id']

gene = pd.read_table('features.old.gz', names = ['gene'])
gene

adata.var_names = gene['gene']

# try do some qc
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
adata

# some stat will be added to obs & vars after qc
adata.obs
adata.var

# add cell meta
cell_meta = pd.read_csv('~/文档/liulab-data-analysis/Archive/covid19/data/ren_zhang2021/GSE158055_cell_annotation.csv.gz')
cell_meta

cell_meta.index = cell_meta['cellName']
cell_meta

new_obs = pd.concat([adata.obs, cell_meta], axis=1)
new_obs

adata.obs = new_obs
adata.obs

# subset B cells
adata_b = adata[adata.obs.majorType == 'B',:]
adata_b

adata_b.write('adata_b.h5ad')

# also save all data
adata.write('ren2021.h5ad')

# %%
adata_b = sc.read_h5ad('data/ren2021/ren2021_b.h5ad',cache=True)

# %%
adata_b.obs
# %%
sample_meta = pd.read_csv('~/文档/liulab-data-analysis/Archive/covid19/data/ren_zhang2021/sample.meta.csv')
# %%
sample_meta
# %%
sample_svr = sample_meta[['title','CoVID.19.severity']].rename(columns={"title":"sampleID"})
# %%
sample_svr
# %%
adata_b.obs = adata_b.obs.merge(sample_svr, how = 'left')
# %%
adata_b.obs.index = adata_b.obs.cellName

# %%
adata_b.var['mt'] = adata_b.var_names.str.startswith('MT-')
sc.pp.calculate_qc_metrics(adata_b, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)

# %%
sc.pl.violin(adata_b, ['pct_counts_mt'],
             jitter=0.4)

# %%
sc.pp.normalize_total(adata_b, target_sum=1e4)
sc.pp.log1p(adata_b)
# %%
sc.pl.violin(adata_b, 'HMCES', groupby='CoVID.19.severity')
# %%
con_deg = sc.tl.rank_genes_groups(adata_b, groupby='CoVID.19.severity',reference='control',method='wilcoxon')

# %%
hmces_deg = sc.tl.rank_genes_groups(adata_b[:,'HMCES'], groupby='CoVID.19.severity',method='wilcoxon')
# %%
hmces_deg.to_df
# %%
